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Particle Swarm Optimization: An efficient method for 
tracing periodic orbits in 3D galactic potentials 



in 
o 
o 

(N 
00 



> 

(N 
O 

in 
o 

S3 

9 L,: 
6 

c3 



X 



Ch. Skokos, 1 ' 2 *! K. E. Parsopoulos, 3 P. A. Patsis 1 and M. N. Vrahatis 3 

1 Research Center for Astronomy, Academy of Athens, Soranou Efesiou 4, GR-11527 Athens, Greece 

2 Department of Mathematics, Division of Applied Analysis and Center for Research and Applications of Nonlinear Systems (CRANS), 
University of Patras, GR-26500 Patras, Greece 

3 Department of Mathematics, University of Patras Artificial Intelligence Research Center (UPAIRC) , 
University of Patras, GR-26110 Patras, Greece 



Accepted; Received; in original form 



ABSTRACT 

We propose the Particle Swarm Optimization (PSO) as an alternative method for 
locating periodic orbits in a three-dimensional (3D) model of barred galaxies. We de- 
velop an appropriate scheme that transforms the problem of finding periodic orbits 
into the problem of detecting global minimizers of a function, which is defined on the 
Poincare Surface of Section (PSS) of the Hamiltonian system. By combining the PSO 
method with deflection techniques, we succeeded in tracing systematically several pe- 
riodic orbits of the system. The method succeeded in tracing the initial conditions of 
periodic orbits in cases where Newton iterative techniques had difficulties. In partic- 
ular, we found families of 2D and 3D periodic orbits associated with the inner 8:1 to 
12:1 resonances, between the radial 4:1 and corotation resonances of our 3D Ferrers 
bar model. The main advantages of the proposed algorithm is its simplicity, its ability 
to work using function values solely, as well as its ability to locate many periodic orbits 
per run at a given Jacobian constant. 
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1 INTRODUCTION 

The periodic orbits of an autonomous Hamiltonian system, 
as well as their stability play a crucial role for the dynamical 
behavior of the system. Orbits that are located near stable 
periodic orbits are ordered, while, near unstable periodic or- 
bits chaotic motion occurs. Therefore, by locating the main 
periodic orbits of a system and following their stability prop- 
erties as one of its parameters changes, we obtain valuable 
information on the ordered or chaotic nature of motion in 
the system. 

The orbital study of barred potentials has pro- 
vided useful i nformation on the structure of galac- 
tic bars fc.f. | Contop oulo£ Il98fj: [Athanassoulal [l983 
Contopoulos fc Grosbe)lll989l ISellwood fc Wilkinson! [l993t 
Pfennigedll996t lPatsisll2004D . [n 2-dimensional (2D) mod- 
els, the galactic bar is supported by regular orbits trapped 
around the so called 'xl' per iodic orbits, which are elongate d 
along the bar major axis (IContopoulos fc Grosb0llll98gD . 
Based on the fact that these orbits do not extend beyond 
the corotation resonance. IContopoulosI il980T) predicted that 
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bars should end at, or before, corotation. In 3D models, the 
planar xl family has in general large unstable parts and, 
thus, its orbits are not sufficient in building the bar. How- 
ever, other families of periodic orbits that bifurcate from xl 
have large stable parts and support the bar. These fami- 
lies build the so-called 'xl-tree' (Skokos, Patsis & Athanas- 
soula 2002a,b). Specific families of the xl-tree are associated 
with certain morphological features observed in real galaxies 
(Patsis, Skokos & Athanassoula 2002, 2003a). Although the 
basic morphological features of barred galaxies are related 
to the presence of orbits of the xl-tree, orbits not belong- 
ing to this tree can also influence the galaxy's morphology. 
Recently, it has been shown (Patsis, Skokos & Athanassoula 
2003b) that inner rings in barred galaxies are associated 
with specific families of periodic orbits influenced by the 
4:1, 6:1 and 8:1 resonances, which are located just beyond 
the end of the bar ( for a definition of the resonances see e.g. 
IContopoulosI J2002D Section 3.1.1). 

The basic families belonging to the xl-tree are usually 
located very easily, but, as we approach corotation, tracing 
periodic orbits that are influenced by high order resonances 
becomes a difficult and challenging problem. The difficulty 
in finding such orbits is mainly due to the fact that at this 
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region there exist many periodic orbits close to each other, 
most of which are unstable. Thus, one is looking for small 
or even tiny islands of stability in a region of the phase 
space of the system, which is mainly characterized by chaotic 
behavior. Usually once a periodic orbit is located, the whole 
family in which it belongs can be found as a parameter of 
the system (e.g., the Jacobi integral) changes. In this way, 
we can follow the morphological evolution and the stability 
transitions of a family and determine its physical importance 
for the system. 

The problem of finding the initial conditions of periodic 
orbits for a given parameter set of a Hamiltonian system is 
the starting point of the present paper. Using an appropri- 
ate scheme we transform the aforementioned problem into a 
minimization problem, where the minimizers of a particular 
function defined on the Poincare Surface of Section (PSS) of 
the system correspond to the periodic orbits of the Hamil- 
tonian. Then, by applying an efficient optimization method, 
Particle Swarm Optimization (PSO), and using deflection 
techniques, we locate periodic orbits, and follow them as a 
parameter of the system varies. This procedure is applied 
successfully on a 3D galactic potential of a Ferrers bar. 

The paper is organized as follows. In Sec. [5] we describe 
the particular model that is used for our orbital calculations. 
Sec^Jis devoted to the thorough description of the proposed 
algorithm. In particular, in Sec. 13. II the procedure of trans- 
forming the problem of detecting periodic orbits into a min- 
imization problem is described, in Sec. l3.2l the PSO method 
for addressing it is presented, while Sec. 13.31 is devoted to 
the deflection technique, which allows us to detect further 
periodic orbits. In Sec. |1J we report the obtained periodic 
orbits and discuss their physical importance. Finally, in Sec. 
|S] we discuss the effectiveness of the proposed numerical 
scheme and present our conclusions. 



2 THE GALACTIC POTENTIAL 

The 3D barred galaxy model that we use is described in 
detail in ISkokos et atl j2002af) . It consists of a Miyamoto 
disk, a Plummer bul ge and a Ferrers bar. The potential of 
the Miyamoto disk iMivamoto fc Nagailll975l) is given by 
the formula, 
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where Md represents the total mass of the disk, A and B are 
scale lengths such that the ratio B /A gives a measure of the 
flatness of the model, and G is the gravitational constant. 
The bulge is a Plummer sphere, i.e., its potential is given 
by, 
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where e s is the bulge scale length and Ms is its total mass. 
Finally, the bar is a triaxial Ferrers bar with density, 



p(m) = 



where, 
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In Eq. Jl}, a, b and c are the principal semi-axes, while Mb 
denotes the mass of the bar component. The corr espondin 
poten tial, Vb , as well as the forces are given in IPfennige: 
( 1984) in a closed form, which is well suited for numerical 
treatment. 

Regarding the Miyamoto disk, we use A = 3 and B — 1, 
and for the axes of the Ferrers bar we set a — 6, b = 1.5 
and c = 0.6. The masses of the three components satisfy 
G(M D + M S + M B ) = 1. In particular, we have GMd = 0.82, 
GM S = 0.08, GM B = 0.10 and e s = 0.4. The length unit is 
taken as 1 kpc, the time unit as 1 Myr and the mass unit as 
2x 1O 11 M0. The bar rotates with a pattern speed, f2(,=0.054, 
around the z-axis, which corresponds to 54 km sec -1 kpc -1 , 
and places corotation at 6.13 kpc. The Hamiltonian govern- 
ing the motion of a test particle can be written in the form, 

H = 7j(pl+p 2 y + pl) + V D + Vs + V B - flb(xp y - yp x ), (5) 

with p x — x — fli,y, p y = y+Qi,x and p z = z being the canon- 
ical momenta. The numerical value of H will be vaguely 
reported as the 'energy', Ej, of the system. 



3 DESCRIPTION OF THE ALGORITHM 

Swarm Intelligence methods are stochastic optimization, 
machine learning and classification procedures that model 
intelligent behavior (Bonabeau, Dorigo & Theraulaz 1999; 
Kennedy & Eberhart 2001). They are closely related to the 
methods of Evolutionary Computation, which consists of al- 
gorithms motivated from biological genetics and natural se- 
lection. A common characteristic of all these algorithms is 
the exploitation of a population of search points that prob e 
the search space simultaneously llSchwefell 1994t lFogell20o'(ih . 

Particle Swarm Optimization (PSO) belongs to the cat- 
egory of Swarm Intelligence methods. The development of 
PSO sprang from the simulation of social dynamics of flock- 
ing organisms, such as insect swarms, which are governed by 
fundamental rules like nearest-neighbor velocity matching. 
In nature, information is communicated among the mem- 
bers of bird flocks and fish schools, enhancing their ability 
to search for food and enabling th em to move synchronized 
without colliding jMillonasI Il994'l . The social behavior of 
animals, and in some cases of humans, is governed by simi- 
lar rules. There is a general belief, and numerous examples 
coming from nature enforce it, that social sharing of infor- 
mation among the individuals of a population, provide an 
evolutionary advantage. 

The dynamics of population in PSO resembles the col- 
lective behavior and self-organization of socially intelligent 
organisms. The individuals of the population exchange in- 
formation and benefit from their discoveries, as well as the 
discoveries of other companions, while exploring promising 
areas of the search space. In our case, the Poincare section is 
our search space, and periodic orbits of a Hamiltonian sys- 
tem are computed through the minimization of a function. 
In the context of function minimization, promising areas of 
the search space are characterized by low function values. 
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3.1 The objective function 

A suitable way to study the stability of orbits in the 6D 
phase space of the Hamiltonian system @ is the well- 
known method of the Poincar e Surface of Section (e.g. 
iLieberman fc Lichtenberel Il992h . Instead of following the 
time evolution of an orbit in the whole phase space, we con- 
fine our study on an appropriately chosen subspace of it. In 
our case, the PSS is the subspace (x,z,x,z) of R 6 , defined 
by the conditions, y — 0, y > 0. The major axis of the bar 
lies along the y axis. Thus, for a given value of the energy 
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(where denotes the transpose of a matrix), on the PSS is 
fully defined as y — and y can be obtained by solving Eq. 
keeping only the positive found value of y. 

In this way, only the 4 initial conditions of an orbit 
on the PSS are necessary for identifying the orbit. Then, 
the time evolution of the orbit is derived by solving the 
Hamilton's equations of motion. The next intersection of 
the orbit with the PSS, y = 0, y > 0, is denoted as, 

*(Xo) = ($4^()),$ z (Xo),$±(X ),$i(Xo)) T :R 4 ^R 4 ,(6) 

which is obviously a point belonging to the 4-dimensional 
PSS of the system. By using the notation O4 = (0, 0, 0, 0) , 
the initial conditions, X, of a p-periodic orbit of the system 
satisfy the equations: 
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Thus, finding the initial conditions, X, of a p-periodic or- 
bit is equivalent to solving Eq. (J7J, which in turn is equiv- 
alent to computing the globa l minimizers of the function 
JParsopoulos fc Vrahati3l2004l . 



f(X) = {§1{X) - xf + m{X) - zf + 

mx) - xf + mx) ~ if . 
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This function is called the objective function, and it is actu- 
ally the square of the Euclidean distance on the PSS between 
the initial point of an orbit and its pth intersection with 
the PSS. Obviously, if the studied orbit is p-periodic this 
distance is zero. Thus, the aforementioned technique trans- 
forms the problem of finding periodic orbits into the problem 
of computing the global minimizers of the function f(X). 

3.2 The PSO method 

Let us now describe the procedure of finding a minimizer 
of the objective function f(X) defined on a 4-dimensional 
search space, SCI 4 , which is a subspace of the PSS. As 
already mentioned, PSO is a population based method, i.e., 
it exploits a population of individuals to probe for promis- 
ing regions of the search space, simultaneously. The popu- 
lation is called a swarm and the individuals (i.e., the search 
points) are called particles. In our case the swarm is a set 
of initial conditions on the PSS and a particle is a point 
X — (x, z, x, z) T on the PSS, which corresponds to an ini- 
tial condition of an orbit of Hamiltonian (JSj- In every it- 
eration of the PSO method, we check whether any of the 



particles is a minimizer of the objective function, f(X), of 
Eq. If it is, the minimizer is recorded, otherwise, each 
particle is moved to a new position in the search space S us- 
ing an adaptable displacement called velocity, retaining also 
a memory of the best position it ever encountered. In our 
case the best positions possess lower function values. 

The particles of the swarm exchange information among 
them. There are two variants of the method with respect to 
the number of particle s that share information: t he global 
and the local variant jKennedv fc Eberhart|l200ll) . In the 
global variant of PSO, the best position ever attained by 
all particles of the swarm is communicated among them. In 
the local variant, each particle is assigned to a neighborhood 
consisting of a prespecified number of particles and the best 
position ever attained by the particles that comprise the 
neighborhood is communicated among them. In the latter 
case, the information attained by the particles is spread in 
the swarm slowly, maintaining high diversity of the parti- 
cles for more iterations of the algorithm than in the global 
variant, which converges faster to the best position. Con- 
sequently, the local variant has better exploration abilities 
while the global variant has better convergence rates (ex- 
ploitation). In the present paper we use the local variant 
of PSO as it is more appropriate for the location of peri- 
odic orbits in the chaotic region near the corotation of the 
bar potential. In this region there are many periodic orbits, 
very close to each other, most of which are unstable. Thus, 
in order to locate such orbits, we need an algorithm that 
performs thorough exploration of the phase space with the 
cost of slightly slower convergence. 

Let us consider a swarm consisting of N particles. Each 
particle is in effect a 4-dimensional vector, 

Xi = {x il ,Xi2,x i3 ,Xi4) T = (xi, Zi, ii, Zi) T e S, 1 = 1,... ,N.(9) 
The velocities of the particles are also 4-dimensional vectors, 



Vi 



(Vil,Vi2,Vi3,Vi4 



,N. 



(10) 



The best previous position encountered by the i— th particle 
is a point in S, denoted by 



Pi = (pu,Pi2,pi3,Pu) e S. 



(11) 



The positions, Xi, and the velocities, Vi, of the particles 
are randomly initialized, following a uniform distribution 
within the search space. The best positions, Pi, are initially 
set equal to Xi. Each particle is evaluated according to the 
objective function f(X) of Eq. JSJ, i.e., the value f{X%) is 
computed for all particles. Obviously, at the initialization 
phase it holds that f(Pi) = f(Xi). 

Let Mi = {Xi- r ,...,Xi-i,Xi,Xi+i,...,Xi + r}, be a 
neighborhood of radius r of the ith particle, Xi (local vari- 
ant). Then, gi is defined as the index of the best particle in 
the neighborhood of Xi, i.e., 



f(P Si )^f(Pj), j = i-r,...,i + r. 



(12) 



The neighborhood's topology is usually cyclic, i.e., the first 
particle Xi is assumed to follow after the last particle, Xn- 
In the general case we face in the present paper, the search 
space is 4-dimensional. In this case we use a swarm con- 
sisting of iV = 20 particles, while the neighborhood of every 
particle has radius, r = 3. For example, the neighborhood of 
the particle X2 consists of the particles, Xig, X20, Xi, X2, 
Xz, X4 and X5. The specific neighborhood size was selected 
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in order to benefit from the exploitation ability of the lo- 
cal PSO variant, while avoiding very slow convergence rates 
implied by smaller neighborhoods. 

If there exist a particle Xj, such that f(Xj) — 0, 
then the initial conditions of a p-periodic orbit are found. 
In the opposite case, we proceed to the next iteration of 
the algorithm, which is to move the particles to new po- 
sitions, taking into account their history so far. This is 
done according to the equations llClerc fc Kennedy! l2002t 
iParsopoulos fc Vrahatisll2002L 120041) 



(9+1) 



X I + Cl n I P, 



3(9) 



X. 



(9) 



+ 



,(9) 



X 



(9) 



X 



(9+1) 



X 



(9) _|_ 



(13) 
(14) 



where i — 1, 2, . . . , N; \ IS a parameter called constriction 
factor; ci and ci are two fixed, positive parameters called 
cognitive and social parameter respectively; n, r2, are ran- 
dom vectors with components uniformly distributed in the 
interval [0, 1]; and q indicates iterations. All vector opera- 
tions are performed componentwise. The objective function 



f(X) is computed again at the new positions, X\ 
the particles, and the best positions, Pi, i — 
updated as follows 
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Then, the new indices, gi, are determined and Eqs. 113H and 
11411 are applied again and so on. We note that in the case 
that the velocity Vi of a particle would move it outside the 
search space S, the particle is actually moved only up to the 
border of S. The algorithm terminates when a user defined 
criterion is achieved, which means that a global minimizer 
is detected. Since in our case the global minimum of the 
objective function f(X) is known a priori to be equal to 
zero, we consider that a minimum of f(X) is located at 
X*, if f(X*) sC 10" 10 . If this criterion is not fulfilled, the 
algorithm stops as soon as a predefined maximum number 
of iterations, g ma x, is reached. In this case, the PSO method 
applied for the particular initial distribution of particles is 
considered to have failed to compute a p-periodic orbit with 
the desired accuracy. In our experiments, g max was set equal 
to 500. 

Let us now discuss the role of the various parameters 
that appear in Eas. 11311 and 11411 . In early versions of PSO 
feberhart fc Kennedvll995T) there was no actual mechanism 
to control the magnitude of the particles' velocities. Thus, 
they could take arbitrarily high values (swarm explosion), 
resulting in divergence of the swarm. For this purpose, a 
positive parameter, V max , was employed as a threshold on 
the absolute value of the velocity's vector components, 

(Vij, if \vij \ < Vmax, 

Vmax ) if Vij <C Vmax j 

Vmax ; if ax . 

for i — 1, . . . , N, j = 1, . . . , n, and n the dimension 
of the search space. In more rec ent versions of PSO 
IClerc fc Kennedvll2002llTreleall2003ft . the constriction fac- 
tor x nas been introduced as a mechanism for constraining 



the magnitude of the velocities. Stability analysis of the algo- 
rithm described by Eqs. I13H and p4t re sults in the followin g 
formula for the determination of v llClerc fc Kennedvl 2002). 



X = 



2k 



|2-0-^2-4^|' 



(15) 



for <f> > 4, where </> = ci + C2, and K = 1. A complete 
theoretic al analysis of the deriva tion of Eq . I|15ll . can be 
found in IClerc fc Kennedy! (120021) and iTreleal J20031) . The 
cognitive parameter, ci, determines the effect of the distance 
between the current position of the particle and its best 
previous position, Pi, on its velocity. On the other hand, 
the social parameter C2 plays a similar role but it concerns 
the best previous position, P gi , attained by any particle in 
the neighborhood. The particular values of the parameters 
that were used in our computations are \ — 0.729 and ci = 
C2 = 2.05. These values are considered optimal default values 
jClerc fc Kennedvll2002l:lTreleall2003h . 

The parameters r\ and 7-2 introduce stochasticity to the 
algorithm. Therefore, if the best positions are very close to 
each other and lie in the basin of attraction of a local mini- 
mizer of the objective function, then r\, hinder the parti- 
cles from getting trapped in that local minimizer, by moving 
stochastically around it. This behavior enables the particles 
to continue searching in potentially better areas of the search 
space, i.e., closer to the global minimizer. Thus, instead of 
moving directly towards the best positions, the particles will 
oscillate around them. The size, N, of the swarm, as well as 
the size of the neighborhood in the local variant of PSO 
can be selected arbitrarily, although it is a common belief in 
evolutionary algorithms that a population size equal from 2 
to 10 times t he dimension of the problem at hand is a good 
initial guess JStorn fc Pricelll997ft . Moreover, the neighbor- 
hood's size shall be problem-dependent. In simple problems, 
larger neighborhoods result in faster convergence without 
loss of the algorithm's efficiency, while, in complicate prob- 
lems with a plethora of local minima, smaller neighborhoods 
are considered a better starting choice. 



3.3 Detecting further periodic orbits through 
deflection 

By applying the PSO method we are able to detect one, in 
general arbitrary, minimizer of the objective function. How- 
ever, in our case, several minimizers of the objective function 
are required, since, in general, there exist many p-periodic 
orbits in the search space, S. Restarting the PSO algorithm 
does not guarantee the detection of a different minimizer. In 
such cases, the deflection technique can be used. This tech- 
nique consists of a transformation of the objective function, 
f(X), once a minimizer, X* , i = 1, . . . , n m i n , has been de- 
tected (Magoulas, Plagianakos & Vrahatis 1997; Parsopou- 
los & Vrahatis 2002, 2004), 



F(X) = J] T^X'.X^fiX), 

i=l 

with, 

Ti(X; X*,A j )=tanh(A i ||X-Jf*||) 



(16) 



(17) 



where Aj, i = 1, . . . , n m i n , are nonnegative relaxation param- 
eters, and n m i n is the number of the detected minimizers. 
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The transformed function has exactly the same minimizers 
with the original f(X), with the exception of X* . Alter- 
native configurations of the parameter A result in different 
shapes of the transformed function. For larger values of A the 
impact of the deflection technique on the objective function 
is relatively mild. On the other hand, using < A < 1 re- 
sults in a function F(X) with considerably larger function 
values in the neighborhood of the deflected minimizer. 

A point to notice is that the deflection technique should 
not be used on its own on a function whose global minimum 
is zero, as is the case for the function f(X) of Eq. ||SJ. The 
reason is that the transformed function, F(X), of Eq. 1161 . 
will also have zero value at the deflected global minimizer, 
since f(X) will be equal to zero at such points. This problem 
can be easily alleviated by considering the function f(X) = 
f(X) + c, where c > is a constant, instead of f(X). The 
function f(X) possesses all the information regarding the 
minimizers of f(X), but the global minimu m is increased 
from zero to c iParsopoulos fc Vrah atis 2004). The value of 
c does not affect the performance of the algorithm and, thus, 
if there is no information regarding the global minimum of 
f(X), it can be selected arbitrarily large. The effect of the 
deflection procedure on the function f(x) = sin 2 (:r) + 0.1, 
at the point x* — ir, is illustrated in Fig.Q In our case, we 
used the parameter values A = 10 4 and c = 0.1 for all the 
computed periodic orbits. 

When the initial conditions of a p-periodic orbit are 
determined, the rest (p — 1) intersections of the orbit with 
the PSS (which can also be considered as initial conditions 
for the same orbit) are obtained through (p— 1) subsequent 
iterations of &{X) of Eq. ©. We determine the stability of 
the p eriodic orbit using e stablish ed techniques (see for ex- 
ample lBrouckejll969l: Hadjidemetrioulll975l: |pfennigeilll984t 
IContopoulos fc MaenenalJll985l:ISkokosll200ll) . The periodic 
orbit can either be stable (S) or unstable. We note that in 
3D Hamiltonian systems we can have three different types 
of instability, namely simple unstable (U), double unstable 
(DU) and complex unstable (A). The repetitive use of the 
deflection technique allows us to find several periodic orbits 
of a certain period, p. The proposed algorithm for the detec- 
tion of periodic orbits is descri bed in pseudocode in Tabled 
JParsopoulos fc Vrahatisll2004l) . 



4 NUMERICAL RESULTS 

Applying the scheme described in Sec. |H] we were able to 
find , apart from the known 1-periodic orbits of Hamiltonian 
© llSkokos et alll2002al : IFatsis et alJl2003bl) . many new 2D 
and 3D orbits of period 1. We also followed these orbits as 
the energy change, registering simultaneously their stability 
transitions and locating their bifurcations. 

To investigate the performance of the algorithm, we 
searched for periodic orbits of multiplicity 1, at the region 
between the radial 4:1 and corotation resonances. In order 
to locate these orbits we used as search space, S, a subspace 
of the PSS, where the values of the initial conditions vary 
in some suitable intervals. Instead of letting all four vari- 
ables, x, z, x, z, to vary simultaneously, which means that 
the search space will be 4-dimensional, we adopted a more 
efficient scheme, described below. 

First we located the planar (2D) orbits of the system on 



the equatorial plain that exist in the 4-dimensional search 
space S. This was done in two phases. In the first phase, we 
located the 2D orbits starting perpendicular to the y — 
axis, having initial conditions of the form (a;, 0,0,0). This 
means that the actual search space in this phase was 1— 
dimensional, as only the x variable of the initial condition 
changed. Confining our search in an 1-dimensional search 
space, also allowed us to decrease the size TV of the swarm 
to TV = 5, using a smaller radius, r = 1, of the neighborhood 
of every particle. All periodic orbits traced in this phase 
were registered and the corresponding initial conditions were 
transformed into maximizers of the objective function f(X) 
through the deflection technique. The search for 2D periodic 
orbits was completed in a second phase by considering or- 
bits starting not perpendicular to the x axis, having initial 
conditions of the form (x,0,x,0). In this stage the actual 
search space was 2-dimensional, while the swarm consisted 
of TV — 10 particles and the neighborhood radius was set 
to r = 2. As the periodic orbits with x = had already 
been found and prevented from been retraced in the previ- 
ous stage of the application of PSO, we located planar orbits 
with x 7^ 0. Again, these periodic orbits were registered and 
prevented from been detected again, through deflection. 

After the detection of the 2D periodic orbits of the 
system on the equatorial plane, we searched for purely 3D 
periodic orbits. Again the search was performed in differ- 
ent phases, keeping the dimensionality of the search space 
as low as possible. Sin ce many 3D periodic orbits found 
bv lSkokos et al.l <l2002aT) have initial conditions of the form 
(x, z, 0, 0) and (x,0,0,z), we first concentrated our tries to 
detect periodic orbits of this form. For such orbits, the actual 
search space of the PSO method is 2-dimensional, since two 
of the four variables are set equal to zero, so we used again 
TV = 10 and r = 2. We note that orbits with z = and i = 
had already been found and, thus, in these phases we de- 
tected actual 3D orbits. Instead of continuing our search by 
letting only one variable be equal to zero and having all the 
possible 3-dimensional search spaces, we preferred to search 
for orbits with initial conditions of the form (a;, z, x, i), be- 
cause the physically most important periodic orbits had al- 
ready been found. Therefore, the complete search for orbits 
with initial conditions of the form (x, z, x, z) was actually 
performed, using TV = 20 and r = 3. 

The scheme of applying PSO in several phases is more 
efficient than the immediate search for orbits with initial 
conditions of the form (a;, z, x, z), because in every phase we 
were confined in search spaces of dimensionality lower than 
4, which is a computationally easier task. Additionally, we 
traced the periodic orbits in a physically meaningful way for 
the particular system, as we found first the planar 2D orbits 
and then the spatial 3D orbits. 

The main periodic orbits that we detected by using the 
PSO method are presente d in the so-c alled 'characteristic' 
diagram (for definition see lContopoulosll2002l . Section 2.4.3) 
shown in Fig. [5] The diagram gives the x coordinate of the 
initial condition of periodic orbits as a function of the en- 
ergy Ej. Since we present in Fig. only planar periodic 
orbits (i.e., z = z = 0), starting perpendicular to the x axis 
(i.e., x — 0) with y — and y > 0, every point defines com- 
pletely the initial condition of the orbit. Orbits on the (x, y) 
plane with i / and orbits not lying on the (x, y) plane 
are not represented on the characteristic diagram of Fig. [5] 
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XXX 

Figure 1. The effect of the deflection procedure on the function f(x) = sin 2 (a;) + 0.1, at the point x* = n, for A = 0.1 (a), A = 1 (b), 
and A = 10 (c). 



Table 1. The proposed algorithm. Note that K denotes the maximum number of the periodic orbits that can be obtained in one run of 
the algorithm. 



Input: 


Hamiltonian H, period p, maximum number of deflections K. 


Step 1 


Set the stopping flag, SF <— 0, and the counter, k <— 0. 


Step 2 


While (SF = 0) Do 




Apply PSO 


Step 3 


If (PSO detected a periodic orbit, X\) Then 




Compute all points X2, ■ ■ • , X p , of the periodic orbit, 




by iterating function <E> (Eq. J(jJ) on Xi. 


Step 4 


If (k < K) Then 




Apply Deflection on Xi, ■ ■ ■ ,X P , and set 




the counter k <— k + 1. 




Else 




Set SF «- 1 




End If 




Else 




Write "No further periodic orbit was detected" 




Set SF <- 1 




End If 




End While 


Step 5 


Report all detected periodic orbits. 



We also note that black curves correspond to stable peri- 
odic orbits, while gray curves indicate unstable periodic or- 
bits without discriminating between their different instabil- 
ity types. In Fig. [5] we plot parts of the characteristic curves 
of the mai n 1-periodic family of the system , the so-called 
xl family iContopoulos fc Papavannopoulosl ll98(Jl . and of 
the known TPatsisetHn* .I l2003bh 2D f and s families, as well 
as the curves of two new families that were found using the 
PSO method. The characteristic curves of these two families 
are very close to each other and to the characteristic curve 
of family s. So, in order to clearly see them we enlarge in 
Fig. |2[ b) the region enclosed in the small rectangular of Fig. 
[3 a). Also, we can see in Fig. |3b) now close to each other 
are the various families of periodic orbits in this region of the 
characteristic diagram. This is the main difficulty that hin- 
ders other classical methods based on Newtonian techniques 
from locating periodic orbits in this region. 

The two new families obtained through the proposed 
approach, exist for energies greater than a minimum en- 
ergy value, similarly to families f and s. Since close to these 
minimum energy values the morphology of these families is 
influenced by the 8:1 and 10:1 resonances, we named these 
families as e and te families, respectively. 



The e family exists for E 3 > -0.19853. In Fig. [3b), 
we see that the characteristic curve of the e family has two 
branches with stable orbits existing at the upper branch. In 
Fig. [3J we plot orbits of the e family for increasing values of 
the energy belonging to the upper branch (orbits (a), (b) and 
(c)), as well as to the lower branch (orbits (d) and (e)) of the 
characteristic. Along the upper branch, we observe a smooth 
transition from a basically 8:1 morphology (orbit (a)) to a 
10:1 morphology (orbit (c)) through an oval-shaped config- 
uration (orbit (b)). We note that, as the energy increases, 
the orbits develop 'corners' along the minor axis of the bar 
(x axis). An analogous transition from a 6:1 to an 8:1 mor- 
phology appears at the upper branch of the characteristic 
of the s family, as can be seen by com paring Figs. joTa). (b ) 
and (c) with Figs. 6(a), (b) and (c) of lPatsis et alJ (12003b'). 
On the other hand, the morphology of the e family at the 
lower branch of its characteristic is always influenced by the 
8:1 resonance (orbits (d) and (e)), and, as energy increases, 
the orbits develop loops (orbit (e)). We also note that the 
unstable orbits of the lower branch of the characteristic have 
almost straight segments parallel to the bar major axis (y 
axis). Similar morphological evolution was observed also for 
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Figure 2. (a) Characteristic diagram of families xl, f , s, e, te. The dash-dotted curve is the section of the zero velocity surface with 
the (Ej,x) plane. In (b) we give the enlargement of the region included in the small rectangular in the upper right side of (a). 



the s family as can be see n by compar i ng Fig s. |5Jd) and (e) 
with Figs. 6(e) and (f) of lPatsis et all d2003bl) . 

A new 2D family is born by bifurcation from the e fam- 
ily at Ej —0.19847. We call it erl following the nomen- 
clature in ISkokos et al ] l|2002alV The erl family has a 9:1 
m orphology (F i g.lH, in analogy to the srl family (Fig. 6(d) 
of lPatsis et al.l <l2003bl) ). which bifurcates from the s family 
and it is influenced by the 7:1 resonance. 

In Fig. El we see the morphological evolution of the 
evl 3D family, which bifurcates from the e family at Ej ~ 
—0.19849. We note that the projections of the evl orbits on 
the (x,y) plane are similar to the planar 2D orbits of the 
e family for the same energies. So the evl orbits appear as 
oval-shaped and evolve to a 10:1 morphology. 

The te family exists for Ej > —0.198438, and it has 
similar evolution to the e family. At the upper branch of 
its characteristic curve, there exist stable representatives of 
the family (orbits (a) and (b) of Fig.[SJ, while, at the lower 
branch, all orbits are unstable (orbits (d) and (e) of Fig. [SJ . 
Family te has a 10:1 morphology (Fig. |HJa) ) , which evolves 
to a morphology influenced by the 12:1 resonance (Fig.lSJc)). 
The orbits belonging to the upper branch of the character- 
istic curve develop 'corners' along the x axis (Fig. EJc)), 
while the orbits of the lower branch have segments parallel 
to the y axis (Figs.|SJd) and (e)). The 2D family terl that 
bifurcates from family te at Ej ~ —0.198434 is influenced 
by the 11:1 resonance (Fig. [7J, while the 3D family tevl, 
which bifurcates from family te at Ej ~ —0.198435, has an 
oval-shaped projection on the (x,y) plane (Fig.|HJ. 



5 DISCUSSION AND CONCLUSIONS 

In this paper, we proposed and applied an algorithm for lo- 
cating periodic orbits in a 3D galactic potential. Our numer- 



ical scheme is based on the transformation of a root-finding 
problem, i. e., solving Eq. (J7J, into a problem of detecting the 
global minimizers of function / (Eq.[SJ . The detection of the 
minimizers is performed by the PSO method. This transfor- 
mation also enables us to use a deflection technique, so that 
many periodic orbits are traced in one run of the algorithm. 
By using the proposed algorithm we were able to trace 2D 
and 3D periodic orbits close to the corotation region in the 
barred galaxy model described in Sec .[5] Our numerical re- 
sults justify the usefulness of the new scheme, as well as its 
advantages, which springs from its ability to locate periodic 
orbits in regions where many periodic orbits coexist very 
close to each other, even when many of them are unstable. 
In addition, even if a periodic orbit is found, Newton-like 
methods cannot guarantee that repeating the search with a 
different initial guess will result in a new periodic orbit. The 
proposed scheme alleviates this problem. 

In order to demonstrate the effectiveness of the pro- 
posed scheme, we discuss now in detail its application on 
the simplest case that we faced in the present paper. This 
case was the location of planar, 2D, orbits of period p = 1, 
starting perpendicular to the y — axis, having initial con- 
ditions of the form (x, 0, 0, 0). Restricting the search for peri- 
odic orbits close to the corotation we set Ej = 0.1984 and let 
x varying in the region 3.5 ^ x ^ 5.5. In this case the search 
space of the PSO method was 1-dimensional and so we used 
a swarm of TV = 5 particles, while the neighborhood of ev- 
ery particle had radius r = 1. Every iteration of the PSO 
method performed the evaluation of the objective function 
f(X) (Eq.|8J with p = 1, for N = 5 times. This means that 
we computed the first intersection with the PSS of TV = 5 
orbits or, in other words, we performed TV = 5 evaluations of 
$(A) (Eq.|SJ, which is the most time consuming computa- 
tion in every iteration of the PSO method. In the particular 
application we asked the tracing of up to 15 periodic orbits 
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Figure 3. Orbits of the 2D e family. From (a) to (b) and then to (c) we see the evolution of the orbital morphology along the upper 
branch of the characteristic of the e family as energy increases, while from (d) to (e) we see the evolution along the lower branch of the 
characteristic plotted in Fig. [5] Orbits (a) and (b) are stable while orbits (c), (d) and (e) are simple unstable. We remind that the bar 
major axis lies along the y axis and the units are in Kpc. 




Figure 4. Orbits of the 2D family erl. From (a) to (b) and then to (c) we see the orbital evolution of the family as the energy increases, 
moving away from its bifurcation point from family e. Orbits (a) and (b) are stable, while orbit (c) is simple unstable. 
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Figure 5. 3D stable orbits of the evl family. For every orbit we plot its projection on the (x,y), (x,z) and [y. z) planes. From (a) to 
(b) and then to (c) the energy increases moving away from the bifurcation point of evl from the e family. 



in one run of the algorithm and we succeeded in locating all 
the planar 2D orbits of families e, te, erl and terl existing 
in the search space, as well as, several banana-like orbits 
(not discussed in the paper) and orbits belonging to the s 
family. All the periodic orbits were located with accuracy of 
at least 1CT 10 . The whole run performed about 1350 itera- 
tions of the PSO method (which means about 90 iterations 
per periodic orbit), involving a total of about 6750 evalu- 
ations of &(X ). A main advantage of our approach is that 
it does not need good initial guesses close to the real initial 
conditions of the periodic orbits it traces. As an example we 
refer to the periodic orbit belonging to the terl family that 
was found by the above mentioned run. This orbit is simple 
unstable and the usual Newton iterative method for finding 
it, converges to its initial condition xo for initial guesses of x 
satisfying \x — xo\ < 7 • 10~ 5 , after about 10 successive iter- 
ations. We note that in the general case of a 4-dimensional 
PSS every iteration of the Newton scheme involves the com- 



putation of the Jacobian matrix of the corresponding set of 
equati ons, which requires 4 eval uations of <&(X) (see for ex- 
ample lPfenniger fc Friedlill993Tl . Thus, in principle, in order 
to be able to locate all the 2D periodic orbits in the interval 
x £ [3.5, 5.5] found by the PSO method, one should compute 
all the orbits with initial conditions on a grid of this interval 
whose grid step should at least be equal to 7 ■ 10 -5 . These 
procedure requires about 28500 evaluations of $(A), not 
taking into account the successive iterations of the Newton 
method to converge to the actual initial condition with the 
desired accuracy. We note that in the case of the terl peri- 
odic orbit this convergence requires about 10 evaluations of 
$(X). So it is evident that the computational effort needed 
by our numerical scheme is significantly lower with respect 
to the Newton iterative method even in the simple case of 
the 1-dimensional search space. 

As we pointed out, the main problem that Newton-like 
techniques face in cases where many periodic orbits coex- 
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Figure 6. Orbits of the 2D te family. From (a) to (b) and then to (c) we see the evolution of the orbital morphology along the upper 
branch of the characteristic of the te family as energy increases, while from (d) to (e) we see the evolution along the lower branch of the 
characteristic depicted in Fig. [2] Orbits (a) and (b) are stable, orbits (c) and (e) are double unstable and orbit (d) is simple unstable. 
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Figure 7. Orbits of the 2D family terl. From (a) to (b) and then to (c) we see the orbital evolution of the family as the energy increases, 
moving away from its bifurcation point from family te. Orbit (a) is stable, while orbits (b) and (c) are simple unstable. 
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Figure 8. A 3D stable orbit of the tevl family. 



ist very close to each other, is that they need a very good 
initial guess of the position of the periodic orbit in order 
to find it. This difficulty of the Newton iterative schemes is 
well known and some efforts to improve their efficiency have 
already been done. As an example, we refer to the paper of 
IPfenniger fc Friedlil <l 19931) where the authors evaluate the 
initial conditions in every successive iteration of the New- 
ton scheme by an appropriate least squares technique. We 
emphasize that our approach is completely different as we 
do not try to find the roots of the system of Eq. J7J but 
we transform this problem into a minimization one. We also 
note that our scheme uses only the values of function $ (Eq. 
|SJ and not its derivatives like Newton iterative techniques 
do. This makes our computation easier and more accurate, 
since the usual method to evaluate the derivatives of func- 
tion $, which is not known in a close analytical form, is to 
approximate them by finite differences. 

Summarizing the main advantages of the proposed algo- 
rithm with respect to Newton-like methods we could men- 
tion that our method is faster, it is simple and can be imple- 
mented easily, it works using function values solely, and also 
it has the ability to locate many periodic orbits per run. 

The orbital behavior of the e and te families and their 
bifurcations, which were detected by t he proposed scheme , 
as well as the behavior of the s family jPatsis et al . 2003b), 
helps us to establish the general behavior of the families ex- 
isting beyond the radial 4:1 resonance gap in our model. The 
morphology of the basic 2D families are influenced by two 
successive even resonances: the s family is influenced by the 
6:1 and 8:1 resonances, the e family is influenced by the 8:1 
and 10:1 resonances and the te family by the 10:1 and 12:1 
resonances. From these families, 2D families influenced by 
the in-between odd resonances bifurcate: srl is influenced 
by the 7:1 resonance, erl by the 9:1 resonance and terl by 
the 11:1 resonance. In addition, the main 3D bifurcations of 
the s, e and te families (svl, evl, tevl families respectively) 
have, in general, projections on the (x, y) plane similar to the 
main families. This is in accordance with the observed fre - 
quency of the various inner rings m orphology jButalll995l) , 
as proposed bv lPatsis et al ] J2003bl) . 
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